Dynamical order, disorder and propagating defects in homogeneous system of 

relaxation oscillators 
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Reaction-diffusion (RD) mechanisms in chemical and biological systems can yield a variety of pat- 
terns that may be functionally important. We show that diffusive coupling through the inactivating 
component in a generic model of coupled relaxation oscillators give rise to a wide range of spatio- 
temporal phenomena. Apart from analytically explaining the genesis of anti-phase synchronization 
and spatially patterned oscillatory death regimes in the model system, we report the existence of 
a chimera state, characterized by spatial co-occurrence of patches with distinct dynamics. We also 
observe propagating phase defects in both one- and two-dimensional media resembling persistent 
structures in cellular automata, whose interactions may be used for computation in RD systems. 
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Understanding how patterns develop in chemical and 
biological systems, e.g., in the course of morphogene- 
sis is one of the enduring challenges of modern sci- 
ence [2|. Investigating pattern formation in systems of 
coupled biochemical oscillators is a promising approach 
towards this objective and has been applied to study, 
for instance, the temporal organization of gene activity 
in cells during the course of development |3(. The pat- 
terns seen in such systems are the dynamical attractors 
that can change with the underlying system parameters 
as the environment evolves with time. Such phenomena 
have been sought to be explained by reaction-diffusion 
models that exhibit patterns such as stripes and spots 
under specific conditions, viz., when an inhibitory chem- 
ical species diffuses faster than an activating species, by 
destabilizing the homogeneous steady state Indeed, 
generalizations of such processes involving differential ex- 
citatory and inhibitory interactions between elements is 
seen in a wide variety of complex systems 0t3|- R- c ~ 
cently such mechanisms have also been proposed as a 
possible basis for computation in biological and chemical 
systems @,@|. 

Time-invariant patterns that are seen in the models 
mentioned above are, however, only a small segment of 
the variety seen in nature, many of which exhibit peri- 
odic activity. Thus, extending the concept of reaction- 
diffusion mechanism to systems of interacting relaxation 
oscillators may allow investigation of spatio-temporal 
patterns in biological systems, where oscillations are ob- 
served across many spatial and temporal scales, e.g., from 
the periodic variations of intracellular molecular concen- 
trations to changes in the activity levels of different 
brain areas ll|. Coherent dynamics of such oscillators 



can produce functionally important collective behavior 
such as synchronization (l2| giving rise to different bio- 
logical rhythms However, synchronized oscillations 
is only one of a number of possible collective phenomena 
that can emerge from such interactions. For example, a 
recent set of experiments on coupled chemical oscillators 



in a microfluidic device [l4| have shown that anti-phase 
synchronization as well as spatially heterogeneous oscil- 
lator death states 15| can occur under different condi- 
tions. Extending the mechanism of coupling by lateral 
inhibition (e.g., via a rapidly diffusing inhibitory chem- 
ical species) to arrays of coupled relaxation oscillators, 
used for modeling biological periodic activity, may reveal 
the underlying mechanism for a variety of spatiotemporal 
phenomena seen in real systems. 

In this paper, we have analyzed in detail a model 
of generic relaxation oscillators (each comprising activa- 
tor and inactivating components) coupled to their near- 
est neighbors through lateral inhibition via diffusion of 
the inactivating component. The model reproduces all 
the spatiotemporal patterns observed in the experiments 
mentioned earlier and its simplicity allows an analytical 
understanding of their genesis. In particular, we have 
given possibly the simplest theoretical demonstration of 
the existence and stability of an anti-phase synchronized 
state for coupled relaxation oscillators. In addition to 
the patterns reported in experiments, we also observe 
other states such as, attractors corresponding to spatially 
co-existing dynamically distinct configurations which we 
term chimera states. Although homogeneous arrays of 
generic relaxation oscillators have been studied exten- 
sively, the observation of such spatially heterogeneous at- 
tractors for these systems is a novel finding of our paper. 
We have performed for the first time a systematic charac- 
terization of the basins of attraction for various patterns 
seen in the model that also demonstrates the unexpected 
robustness of the chimera states and suggest that all the 
observed states can be obtained in suitable experiments. 
In addition, we report the occurrence of phase defect- 
like discontinuities moving ballistically through the sys- 
tem that on collision with each other can produce com- 
plex patterns. We observe analogous structures in two- 
dimensional media that have a striking resemblance to 
persistent configurations in cellular automata (CA), e.g., 
"gliders" in the "Game of Life" CA [3] , which have been 



2 



Passive Oscillators 
# — • — • — • — • 

12 3 4 

b c 





Passive 
■• — • — • — • — # 

N-3 N-2 N-1 N N+1 

d e 



3 1 



CD 

E 3 




0.6 




FIG. 1: (color online). Spatio-temporal evolution of a 1- 
dimensional array of coupled relaxation oscillators (iV = 10) 
with passive elements at the boundaries [model system shown 
schematically in (a)] . Pseudocolor plots of the activation vari- 
able u indicate different regimes characterized by (b) syn- 
chronized oscillations (SO), (c) anti-phase synchronization 
(APS), (d) spatially patterned oscillation death (SPOD) and 
(e) chimera state (CS), i.e., co-occurrence of spatial patches 
with dynamically distinct behavior. 



linked to the universal computation capabilities of such 
systems 17|, [l8[ . The observation of these patterns in our 
model system is remarkable considering the simplicity of 
the underlying dynamics. 

The system we have analyzed comprises N relaxation 
oscillators interacting with each other in a specific topol- 
ogy. For the dynamics of individual relaxation oscillators 
we have used the phenomenological FitzHugh-Nagumo 
(FHN) equations, which is a generic model for such sys- 
tems. Each oscillator is described by a fast activation 
variable u and a slow inactivation variable v. 



u = /(it, v) = u (1 — u) (u — a) — v, 
v = g(u, v) = e (k u — v — b), 



(1) 



where, a = 0.139, k = 0.6 are parameters describing 
the kinetics, e = 0.001 characterizes the recovery rate 
of the medium and b is a measure of the asymmetry of 
the oscillator (measured by the ratio of the time spent 
by the oscillator at high and low value branches of u). 
Parameter values have been chosen such that the system 
is in oscillatory regime. Small variations in the values 
do not affect our results qualitatively. To investigate 
spatial patterns generated by interaction between the 
oscillators, they are arranged in a 1-dimensional chain 
[Fig. Q] (a)]. In actual chemical experiments, the beads 
containing the reactive solution are suspended in a chem- 
ically inert medium which allows passage of only the in- 
hibitory chemical species [l4[- In our model, the oscilla- 
tors are diffusively coupled via the inactivation variable 



v. The boundary conditions for the chain take into ac- 
count the inert medium by including non-reactive pas- 
sive elements at each end that are diffusively coupled to 
the neighboring oscillators. The inert medium between 
the oscillators are not considered explicitly, their volume 
being relatively small compared to the reservoirs at the 
boundary. We have verified that inclusion of interme- 
diate non-reactive cells diffusively coupling each pair of 
oscillators do not affect the fixed-point equilibria of the 
system or their stability, subject to suitable rescaling of 
the diffusion constant. The dynamics of the resulting 
system is described by 



v = D v (vi - wo), 

Ui = f(Ui,Vi), 

ii = g(ui, Vi) + D v (i>i_i + v i+ i - 2 Vi), 



(2) 
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where i — 1,2, ... ,N and the diffusion constant D v rep- 
resents the strength of coupling between neighboring re- 
laxation oscillators through their inactivation variables. 
For most results reported here N — 10 although we have 
used higher values of N upto 1000 to verify that the qual- 
itative results are not sensitively dependent on system 
size. We have also verified that the boundary conditions 
do not sensitively affect the results by considering pe- 
riodic boundaries and observing patterns qualitatively 
identical to those reported here. The dynamical equa- 
tions are solved using an adaptive Runge-Kutta scheme 
(e.g., rkf45 19]). The behavior of the system for each set 
of parameter values b and D v is analyzed over many (10 4 ) 
initial conditions, with each oscillator having a random 
phase chosen from a uniform distribution. 

Fig. [TJ (b-e) shows a variety of asymptotic spatio- 
temporal patterns that are observed in the model sys- 
tem: synchronized oscillations (SO) with all elements 
(except those at the boundary) having the same phase 
(b), anti-phase synchronization (APS) with neighbor- 
ing elements in opposite phase (c), Spatially Patterned 
Oscillator Death (SPOD) regime where the oscillators 
are arrested in various time- invariant states (d) and 
Chimera States (CS) where oscillating regions co-exist 
with patches showing negligible temporal variation (e). 
However, these do not exhaust the range of possible 
spatio-temporal phenomena that are observed including 
propagating structures that are discussed later. Note 
that, both APS and SPOD states have been observed 
experimentally in chemical systems. Although the latter 
have sometimes been referred to as "Turing patterns" 
in the literature, SPOD is distinct as it is not obtained 
through destabilization of a homogeneous equilibrium 
(Turing instability) but through a process of oscillator 
death. 

To investigate the robustness of the observed patterns 
in detail, we numerically estimate the size of their basins 
of attraction in the (b, D v ) parameter space (Fig. [5]). This 
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provides an indication as to whether a pattern can be 
observed in actual experiments as for this they must be 
immune to the unavoidable noise present in any experi- 
ment. In order to segregate the distinct pattern regimes 
in the (b, D v ) space [Fig.[2](a)] we introduce the following 
order parameters. The number of non-oscillating cells in 
the bulk of the system, N no , i.e., cells for which the vari- 
ance of the activation variable u, of (uj), is zero, is used 
to characterize the SPOD (N no = N) and CS regimes 
(0 < N no < N). The SO and APS states both have 
all elements in the bulk oscillating. However, SO is dis- 
tinguished by having all oscillators in the same phase as 
measured by the variance of the activation variables u, 
{af(u)) t = 0, where ( ) t represents time average. We can 
also define the synchronization among the oscillators in 
two distinct sub-lattices, namely, those which occur at 
even-numbered and those which occur at odd-numbered 
sites of the chain, as measured by the time-averaged 
variance of the activation variable, viz., (c^ ven {u))t and 
( a odd( u ))t- This pair of order parameters are zero for 
both SO and APS states; however, if (of (u)) t > 0, it sig- 
nifies the APS regime. In practice, the different regimes 
are characterized by thresholds whose specific values do 
not affect the qualitative nature of the results. Fig. [5] (a) 
indicates the parameter regions where SO, APS, SPOD 
and CS states are observed for more than 50% of initial 
conditions (i.e., they have the largest basin). As already 
mentioned, the system also exhibits other regimes apart 
from the above ones, which occur in regions of (6, D v ) 
parameter space shown in white. 

While diffusive coupling in a homogeneous system of 
oscillators is expected to promote the SO state |2(|, a 
striking result from this phase diagram is that for cer- 
tain parameter values the APS state has a very large 
basin of attraction [Fig. Hth)]. The existence of APS 
is somewhat counter-intuitive as for diffusively coupled 
identical isochronous oscillators the only stable attrac- 
tors are synchronized oscillations or oscillator death [2(| ■ 
While anti-phase synchronization has been seen earlier in 
a pair of identical oscillators 2l| , it is not obvious that 
APS will have a large basin of attraction for an array of 
oscillators. To understand the origin of such anti-phase 
oscillations we consider a simple model that captures the 
essence of relaxation oscillation phenomena and can be 
solved exactly. We consider the relaxation limit (e-)0 
in FHN system) and extreme asymmetry where the limit 
cycle has a slow segment in which the system spends the 
entire duration of the oscillation period (the remaining 
segment of the cycle being traversed extremely fast). Un- 
der these considerations, we obtain the one-dimensional 
dynamical system: x = iv(x), where x parameterizes the 
slow part of the limit cycle and can be redefined to be- 
long to the interval (0, 1). Fig. [3] (a) shows a schematic 
diagram of the trajectory of the limit cycle, where the 
system spends almost its entire oscillation period on the 
solid branch (the return from x = 1 to x = 0, shown 
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FIG. 2: (color online). Different dynamical regimes of a 1- 
dimensional array of coupled relaxation oscillators (N — 10) 
in the D v — b parameter plane showing regions where the ma- 
jority (> 50%) of initial conditions result in synchronized os- 
cillations (SO), anti-phase synchronization (APS), spatially 
patterned oscillator death (SPOD) and chimera state (CS). 
(b) Variation of the attraction basin size for the different 
regimes mentioned above (measured as fraction of initial 
states reaching the attractor) with coupling strength D for 
b = 0.064 [i.e., along the broken line shown in (a)]. In practice, 
the regimes are distinguished by thresholds applied on the or- 
der parameters of («i)> (of(u)) t , (of«en(«))t and (<r 2 odd (u)) t , 
which have been taken to be 0.05 for the present figure. Basin 
sizes have been estimated using 10 4 initial conditions. 



by the broken line, is considered to be instantaneous). 
The model can be exactly solved if cj(x) is a constant 
(= Li, say), although the geometrical argument is valid 
for any arbitrary positive definite function defined over 
the interval (0,1). By appropriate choice of time scale, 
we can set the period uj^ 1 = 1 without loss of generality. 
A system of two such diffusively coupled oscillators can 
be described by 



±i = 1 + D (x2 — xi), ±2 = 1 + D [x\ - x 2 ) 



(3) 



Given the values of x\,X2 at some arbitrary initial time 
t', the solution of Eqn. ((3J at a later time t follow the 
relations: 



xi(t) 
xi(t) 



x 2 (t) 
x 2 (t) 



Xx{t') + X2{t') - 
[ Xl (t') - X 2 (t')} 



■2(i- 
exp[ 



■0, 

-2D(t 



(4) 



till time t" when max(xi, X2) reaches x = 1. After this 
the larger of Xi,X2 is mapped back to x = (because of 
the instantaneous nature of the remaining segment of the 
limit cycle) and t' in Eqn. (|4|) is replaced by t" . Using 
the above exact solution of the coupled system (J3j) , its 
Poincare map P(x) is constructed in two steps. First, if 
xi starts at and x 2 starts at some point < x < 1, 
we find the location of Xx[= f{x)] at some time t when 
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FIG. 3: (color online), (a) Schematic diagram of a limit cy- 
cle trajectory for an oscillator in the relaxation limit (e — > 0) 
and extreme asymmetry (for details see text) such that the 
oscillator is on the solid line (0 < x < 1) for its entire pe- 
riod, (b) Time-series of two such coupled oscillators [Eqn. Q 
with D = 1] and (c) the Poincare map for the system at 
different coupling strengths D showing stable anti-phase syn- 
chronization, (d) Phase-plane diagram indicating the general 
mechanism (see text) for oscillator death in a system of two 
coupled oscillators (1 and 2). 



X2 = 1 (which is then immediately mapped to X2 = 0). 
Now, starting with X2 — and x\ = f{x), when x\ = 1 
we find the location of X2- x' = f(f(x)) — P(x). Using 
Eqn. gj, with Xi[t') = 0, x 2 {t') = X, x^t) = f(x) and 
x 2 (t) = 1, we get 



f{x) = 1 + D- 1 W[-Dx exp{D(x - 2)}], 



(5) 



where W is the Lambert W-function. Fig. [3J (b) shows 
the Poincare map P(x) = f(f(x)) for different values of 
the coupling strength D. The map has one stable and one 
unstable fixed point, which correspond to the anti-phase 
synchronized (APS) and synchronized oscillating (SO) 
states, respectively. Thus, for the model ([3]) we find that 
APS state not only exists but is the only stable state. 
Relaxing the extremal conditions under which this was 
derived may allow a stable SO state to coexist with the 
stable APS state (22J. The derivation we have shown here 
is possibly the simplest one exposing the fundamental 
mechanism for generating APS states in any system of 
diffusively coupled oscillators that can exhibit anti-phase 
oscillations. 

When the coupling D v between oscillators in the ar- 
ray is increased to very high values, we observe that the 
oscillatory regimes (e.g., SO and APS) are replaced by 
time-invariant spatial patterns such as SPOD (Fig. [2|. 
To understand the genesis of SPOD at strong coupling, 
we can again focus on a pair of coupled relaxation os- 
cillators in the relaxation limit (e — > 0). The parameter 
b is chosen such that the u-nullcline is placed symmet- 



rically between the two branches of the u-nullcline with 
the oscillator taking equal time to traverse each branch 
[Fig. [3] (d)]. When the two oscillators (1 and 2) are in 
opposite branches (as shown in the schematic diagram), 
the two opposing forces acting on each oscillator, corre- 
sponding to the coupling [Fd = D v (v2 — v\)] and the in- 
trinsic kinetics (F n ) respectively, can exactly cancel when 
the coupling is strong resulting in oscillator death. Note 
that only the component of Fd parallel to the intrinsic 
force F n (F d ) needs to balance F n as F^ has no effect 
in the relaxation limit. The symmetry of the oscilla- 
tor ensures that the force due to the intrinsic kinetics 
(F n ) for the two oscillators are identical but oppositely 
directed in the steady state. The occurrence and stabi- 
lization of this heterogeneous time-invariant state (as 1 
and 2 need to be in different branches, corresponding to 
very different values of «i,2) is the key to the occurrence 
of SPOD at strong coupling. At intermediate values of 
coupling D v between oscillators in a large array, the com- 
petition of this mechanism with the natural oscillatory 
dynamics (seen at low coupling) may give rise to chimera 
states that exhibit characteristics of both, i.e., contain- 
ing patches of elements that are oscillating as well as 
segments that show SPOD-like pattern. This CS regime 
is especially interesting as the system exhibits a strik- 
ingly heterogeneous dynamical state in spite of the bulk 
of the array being homogeneous with identical elements 
coupled in the same fashion. We have explicitly verified 
that the occurrence of CS is not dependent on bound- 
ary conditions by reproducing it also in systems with 
periodic boundaries. The observation of such states in a 
generic model of relaxation oscillators suggests that they 
should be present in a wide class of systems, and indeed 
similar phenomena have been recently reported in a spe- 
cific chemical system model (23|. Note that, the chimera 
state described here comprises regions with dynamically 
distinct behavior, and is different from its namesake that 
refers to co-occurrence of coherent and non-coherent do- 
mains [2^ | . 

Apart from the spatio-temporal patterns shown in 
Fig. Q] (b-e) we also observe attractors corresponding to 
point-like "phase defects" (i.e., there is a discontinuity 
of phase along the oscillator array at this point) moving 
in the background of system-wide oscillations. As seen 
from a typical example of such states [Fig. H] (a)], after 
initial transients these defects move in the medium with 
interactions between two such entities resulting in either 
the two getting deflected in opposite directions, or either 
one or both getting annihilated (unlike defects in non- 
oscillatory media such as domain walls which annihilate 
on collision 
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While the boundary for systems with 
passive elements at the ends is a source of new defects 
entering the medium, similar persistent structures are 
also seen in systems with periodic boundary conditions 
where, in the steady state, a conserved number of defects 
can reflect off each other indefinitely [Fig. @] (b)]. 
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FIG. 4: (color online), (a-b) Spatio-temporal evolution of 
a system of coupled relaxation oscillators showing traveling 
waves of phase defects in (a) a linear array with passive ele- 
ments at the boundaries and (b) with periodic boundary con- 
ditions, (c-d) Propagating defects in two-dimensional media 
with periodic boundary condition showing (c) two horizon- 
tally moving "gliders" and (d) collision of two "gliders" . For 
clear visualization of the motion of the spatially extended de- 
fects, snapshots of the two-dimensional medium are taken at 
intervals which are multiples of the oscillation period for the 
mean activity of the system r. Animations are available in 
Ref. [13]. 



To observe how these propagating defects manifest in 
higher dimensional systems, we consider a 2-dimensional 
array of coupled oscillators defined on a torus [Fig. 0] (c- 
d)]. The system can have extremely complicated tran- 
sient phenomena, and for simplicity we show here only 
its asymptotic behavior. For a square lattice, we observe 
that there is a specific configuration of four sites that 
persists indefinitely (reminiscent of the glider configura- 
tions in the 2-dimensional CA "Game of Life" [l6j) and 
can move in horizontal or vertical directions. Interaction 
of such "gliders" can produce complex spatio-temporal 
patterns. Fig. [4] (d) shows two "gliders" that on colli- 
sion move off in directions perpendicular to their incident 
ones (iij ]. 

To conclude, we have shown that a simple model of 
relaxation oscillators coupled via lateral inhibition can 
exhibit a wide variety of striking spatio-temporal pat- 
terns. Our comprehensive investigation has revealed the 
global features of the dynamics and robustness of the 
attractors. As our results are based on a very generic 
model, it suggests that the patterns may be observed in 
a range of experimental realizations, such as electronic 
circuits implementing relaxation oscillators 27| and Pt 
wire undergoing CO oxidation where the system is in 
an oscillatory regime [Hj], apart from microfluidic chem- 
ical systems mentioned earlier. Our initial exploration of 
propagating configurations in 2-dimensional media sug- 



gests that systems of higher dimensions may yield even 
more striking discoveries. It is intriguing to explore the 
possibility of using the propagating defects for computa- 
tion [1, E| , as analogous entities have been used to con- 
struct logic gates in CA [16J. This may tie together two 
of the groundbreaking ideas of Alan Turing, viz., pattern 
formation through reaction-diffusion mechanism and uni- 
versal computation [29|. 

We thank Bulbul Chakraborty and Sudeshna Sinha for 
helpful discussions and the HPC facility at IMSc for pro- 
viding computational resources. 
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